sigma1 = [2,1; 1,2];
sigma2 = [1,0; 0,1];
sigma3 = [1.5,0; 0,1.5];
mu1 = [0;3];
mu2 = [0;0];
mu3 = [5;0];

X1 = gauss_sample(mu1, sigma1, 100)';
X2 = gauss_sample(mu2, sigma2, 100)';
X3 = gauss_sample(mu3, sigma3, 100)';
Xall = [X1; X2; X3];
yall = [ones(100,1); 2 * ones(100,1); 3 * ones(100,1)];

num_components = 2;
gmm = gmm_train(Xall, yall, num_components);
yall2 = gmm_classify(Xall, gmm);
ycorrect = yall == yall2;

figure
hold on
scatter(X1(:,1),X1(:,2),'x','b');
scatter(X2(:,1),X2(:,2),'x','k');
scatter(X3(:,1),X3(:,2),'x','c');
plot_gauss_2d(mu1, sigma1, 1, 'b');
plot_gauss_2d(mu2, sigma2, 1, 'k');
plot_gauss_2d(mu3, sigma3, 1, 'c');

figure
hold on
scatter(Xall(ycorrect,1),Xall(ycorrect,2),'x','g');
scatter(Xall(ycorrect == 0,1),Xall(ycorrect == 0,2),'x','r');
plot_gauss_2d(mu1, sigma1, 1, 'b');
plot_gauss_2d(mu2, sigma2, 1, 'k');
plot_gauss_2d(mu3, sigma3, 1, 'c');

figure
hold on
scatter(Xall(ycorrect,1),Xall(ycorrect,2),'x','g');
scatter(Xall(ycorrect == 0,1),Xall(ycorrect == 0,2),'x','r');
for k = 1 : num_components
    plot_gauss_2d(gmm.model{1}.mu{k}, gmm.model{1}.sigma{k}, 1, 'b');
    plot_gauss_2d(gmm.model{2}.mu{k}, gmm.model{2}.sigma{k}, 1, 'k');
    plot_gauss_2d(gmm.model{3}.mu{k}, gmm.model{3}.sigma{k}, 1, 'c');
end
